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Abstract 

The behavior of weakly coupled self-sustained oscillators can often be well described by phase equa- 
tions. Here we use the paradigm of Kuramoto phase oscillators which are coupled in a network 
to calculate first and second order corrections to the frequency of the fully synchronized state for 
nonidentical oscillators. The topology of the underlying coupling network is reflected in the eigen- 
values and eigenvectors of the network Laplacian which influence the synchronization frequency 
in a particular way. They characterize the importance of nodes in a network and the relations 
between them. Expected values for the synchronization frequency are obtained for oscillators with 
quenched random frequencies on a class of scale-free random networks and for a Erdos-Renyi ran- 
dom network. We briefly discuss an application of the perturbation theory in the second order to 
network structural analysis. 

PACS numbers: 05.45.Xt, 64.60.aq 
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I. Introduction 



The collective behavior of ensembles of interacting units is one of the main topics in complex 
system theory. Different parts of a complex system can be identified as subsystems and 
studied individually, while the interaction between these can lead to emergent properties of 
the whole system. In particular synchronization, the adjustment of internal timescales in 



oscillatory systems which interact locally or through a complex network 
in biological and technical applications 



compl 
0,0, 
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2| , is ubiquitous 



12l |. Recently also chemical 



reactions with feedback control have been proposed to realize specific interaction topologies 
0, Q] • Synchronization can orchestrate macroscopic spatio-temporal periodicity even if the 
individual units are very different from each other and a simple linear superposition of their 
output would be incoherent. While this is a desirable effect in many applications, such as 
coupled Josephson junctions or laser arrays 0, 01 it can also lead to pathological states lik e 



epilepsy or Parkinson's disease [4'] or it can be disastrous when it occurs in constructions 



15| 



The onset of synchronization for very heterogeneous systems has been described as a 
second order phase transition in the limit of large system sizes Above a critical 

coupling strength or below a critical heterogeneity the incoherent state becomes unstable 
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19|. For identical, possibly chaotic. 



20 



2l|. It is known that the 



and global collective behavior can be observed 
subsystems complete synchronization can be possible 
spectral properties of the coupling network play an important role in the transition to 
synchronization and the stability of complete synchronization {20! . I21I . But many 

studies on synchronization in networks have mainly been concerned with the estimation of 
a few important eigenvalues of the network Laplacian 22 1. 



In this paper we study the synchronization frequency in networks of weakly noniden- 
tical, autonomous oscillators with attractive coupling. Under these conditions the 
Kuramoto phase equations (KPE) [17i] can be used to describe the system qualitatively and 
quantitatively. The KPE show a rich collective behavior with transitions from complete 
desynchronization, where the phases are uniformly distributed, to partial synchronization 



with a unimodal distribution of phases or even clustering 



13, 



18|, 



19|, 



24 



25l ] and finally 



frequency synchronization or phase locking, where the phase difference between any two 



2 



oscillators is constant. 



In systems of identical phase oscillators with attractive coupling complete synchro- 
nization is a stable solution of the KPE. We will quantify the frequency heterogeneity of 
the oscillators and derive a perturbation expansion around the well known synchronization 
manifold for identical oscillators in powers of the heterogeneity. In analogy to perturbation 



theory for the continuous, nonlinear Kuramoto Phase Diffusion Equation [26|, in Sections 
II and III we will show two approaches which lead to the same first and second order per- 
turbation terms. In random networks, the expected second order perturbation term of the 
synchronization frequency can be interpreted as a mean value with respect to the spectral 
density of the network Laplacian. Using a random network model for which the spectral 
density of the Laplacian is known, we explicitly calculate the expected synchronization 
frequency in Section IV. We verify out theory by numerical simulations. 

The Kuramoto model 

et us briefly review the Kuramoto Phase Equations (KPE) for discretely coupled oscillators 
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16|, [iTj The dynamics of an ensemble of N autonomous oscillators may be given as 



N 

nm (Xm, X„) , (1) 

m=l 

where X„ defines the state of the oscillator labeled with n, the velocity field F„ allows 
for stable limit cycle oscillations and V^nm describes the coupling between two oscillators 



depending on their state. In his monograph 



17( 1 Kuramoto considered the heterogeneity in 



the oscillators as well as the coupling as a perturbation of a common oscillator dynamics 
F„ = F+5F„. For this common dynamics one can define a uniformly evolving phase variable 
in a neighborhood of the limit cycle. The dynamics of the phases 0„ in linear response to 
the perturbation corresponds to the phase model introduced by Winfree 

N 

(2) 



UJ 

m=l 



Here uj is the natural frequency and Z is called the phase response function of the common 
oscillator dynamics. If the phase differences change slowly over the time of one oscillation 
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then one can use phase averaging techniques [l|, ll7|] to obtain effective phases and phase 
equations which only depend on the phase differences. If we finally assume that the func- 
tional form of the coupling between any two oscillators n and m only differs in a coupling 

constant Anm we obtain the KPEs 

N 

i}n = ^n + ^^ Anm g{'&m " "^n) ■ (3) 
m=l 

The phase coupling function g{Ai)) is periodic. For diffusive coupling it vanishes at zero. 
We assume a positive derivative at zero and approximate the coupling function by its lowest 
Fourier modes as 



g{A^) = sin + 7 (1 - cos A^) . (4) 

The parameter 7 breaks the symmetry of the phase coupling function and can directly be 
associated with the amplitude dependence of the phase velocity in complex Stuart-Landau 
oscillators i.e. a third-order nonlinear effect in the normal form of a supercritical Hopf 
bifurcation also known as nonisochronicity. The effect of nonisochronicity on the ability of 
a system to synchronize and on the formation of sp atio-temporal patterns has been noted 



early on [27] and again stressed recently 



of analytic simplicity 



m 
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2^ whereas it is often disregarded in favor 



18|. 



A fully phase locked state is reached when the oscillators can arrange their phases in 
a way that due to an exact balance of nonidentical natural frequencies and coupling forces 
all oscillators have the same synchronization frequency 

N 

^ = ar]n + '^ Anm g{'&m " ^?n) • (5) 
m=l 

The frequencies ?7„ in this equation are normalized to have unit variance. Then the hetero- 
geneity of the oscillators is quantified by the variance var(co') = of the natural frequencies 
in the system. The mean frequency Co does not necessarily depend on the heterogeneity a 
but here we choose a co-rotating frame of reference where uj = af]. For identical oscillators 
(cr = 0) complete synchronization with identical phases i!}n^ = "^m^ for all n and m, and 
synchronization frequency fi*^^-' = is a solution of Eq.Q with 

N 

= = A_ ^(^(^) - (6) 

771=1 
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Stability of the synchronized state 



Under some weak conditions on the coupling topology one can show that the state of 
complete synchronization is stable. But it has been shown recently, that in a heterogeneous 
coupling network and for large nonisochronicity 7 the stable state of complete synchroniza- 
tion can co-exist with a dynamical equilibrium of complete desynchronization or partial 
synchronization [23]. Conversely, if the nonisochronicity is not too high and the network is 
well connected, complete synchronization is the typical result from random initial conditions. 



A sufficient condition for the stability of complete synchronization of identical oscil- 
lators is that all values Anm are non-negative and the corresponding weighted network is 
strongly connected, i.e. there exists a path between any two nodes. To see this, one can 
consider small deviations (fn from the synchronized solution. Linearizing Eq. ([3]) for a = 
and small deviations around i^^^-* one obtains 

N N 

0n = ^ Anm {^m - ^n) = ^ Lnm V^m , (7) 
m=l m=l 

with the network Laplacian L defined as 

N 

L'nm -Afim ^nm ^ ^ • (8) 

1=1 

Since all row sums J2m ^nm are zero at least one eigenvalue Aq of the network Laplacian 
is also zero, corresponding to a constant shift of all phases along the synchronization 
manifold. If all values Anm are non-negative then it follows from the Gershgorin cir- 
cle theorem that the network Laplacian has only non-positive eigenvalue real parts 
= Aq > ReAi > ■ ■ • > ReA7v_i, where is the number of oscillators. Complete 
synchronization is unique up to a global phase shift, only if the second largest eigenvalue 
real part ReAi is strictly smaller than zero. 



Associated with the relaxation to synchronization is a diffusion process in the oppo- 
site direction of the coupling. If all off-diagonal elements are non-negative, the transposed 
Laplacian L^" can be viewed as a matrix of transition rates for a master equation P = L'^'P 
with a probability vector P. The eigenvalue Aq is non-degenerate if and only if the 
stationary probability distribution Pq is unique, i.e. independent of the initial condition. 

5 



Note, that a strongly connected network of transition rates is sufficient but not necessary 
for that [30,]. In the following we will assume that ReAfc < for all A; > 0. 

II. Perturbation Approach 1 

The algebraic equation Eq.([5]) implicitly defines the synchronization frequency and the 
phases in synchronization (up to global phase shift), even for non-zero heterogeneity. How- 
ever, a stable phase locked solution of Eq.(l5]) or any solution at all may not exist. Only 
for small heterogeneity we can expect that a stable solution exists, that it is close to the 
solution for identical oscillators (a = 0) and that it can be expanded in powers of o as 

(9) 



In this section we follow closely the procedure outlined in [26(] to derive the perturbation 
expansion of the Kuramoto Phase Equations in synchronization Eq.(l5]). We directly insert 
the ansatz Eq.Q into Eq.^S]), use a Taylor expansion of the coupling function around zero 
and regroup the terms according to powers of a. This procedure requires sorting of infinite 
summations and some careful consideration of the index limits. It is shown in detail in the 
Appendix A. However, the result takes a simple form in vector notation 

n«i = (L^(') + b«) . (10) 

Here fl^^^ is the l^^ order perturbation term of Q in Eq.([9]), the vector i?*-'^ is the corre- 
sponding perturbation term for the phases in synchronization, 1 is a constant vector with 
unity entries, the matrix L is the Laplacian of the network, as defined in Equation Eq.([8]) 
and b^') is a vector which depends nonlinearly on all perturbation terms of order lower 
than / (see Eqs. (1141) - (|T6l) ). Equation ( fTOl) can thus be solved iteratively for each order of 
perturbation. In practice, while the amplitude of the terms b^'^ is as small as 0((t'), the 
analytic expression and the expense for its calculation blows up quickly. 

Let us consider a complete, orthonormal set of left and right eigenvectors and pfc 
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of the network Laplacian with 

(11) 

Af-l 



k=0 



and in particular 



Po = 1 and it Po = 1 . (12) 

The left eigenvector Pq, which is the stationary ^solution of the master equation P = Ltp, 
assigns a weight to each node of the network [31]. The solution of Equation f[TI]|) is 

fi« = pt b« , 

(13) 

^ = -Pk- 

Using the short notations Qq = g"{0), Qq = g"'{0) and -^mn = 'dm —"dn^ the first three vectors 
b(i), b(2) and b(3) are 

h^n^ = Vn , (14) 
^ 1 

b^:^ = E ^nm-g'Ml" , (15) 
m=l 

b^:^ = E (aXi^^l + Ig'^'^^r^j) ■ (16) 

m=l ^ ^ 

Equations (|T3 l) -( |T6l) give the first three perturbation terms of the synchronization frequency 
and the relative phases in synchronization. In the next section, we will derive the first 
and the second order terms again, but in a slightly different form which allows for a better 
analysis. 

III. Perturbation Approach 2 

The second order correction Eqs. (fT3|) and (fTS!) of the synchronization frequency depends 
on the second derivative = g"{0) of the phase coupling function at zero. If we are only 



interested in the first and second order perturbation terms we have the freedom to choose 
a different couphng function g{(p) in the equation Eq.Q with the same first and second 
derivative at zero as g{(f) which may be more suitable for an analysis. For the continuous 
Kuramoto Phase Diffusion equations it is known that a non-hnear Cole-Hopf transformation 
= 7~^logp changes the equations in synchronization into an eigenvalue problem of a 



stationary, linear Schrodinger equation 17|, |26|, |27|, |28| . With the same procedure in mind 



we will define an auxiliary coupling function g{(p) as 



m = - (e^'^ - 1) = g{^) + 0{ip') . (17) 
7 

After the transformation 

■dn = - logPn , (18) 

7 

we can bring equation Eq.Q with g{'^) as coupling function into the form of an eigenvalue 
problem 

- EoP = l^V = [ic^^v + L] P = - Hp , (19) 

where the vector p has the entries pn, V,, = diag(T7) is the diagonal matrix of frequencies 
and L is the network Laplacian (Eq.®). This equation has the form of a stationary discrete 
Schrodinger equation for the ground state of a particle hopping between the vertices of the 
coupling graph with the on-site potentials —rin and ground state energy Eq = —'-/Q. If the 
coupling network is symmetric the Hamiltonian H is symmetric, as well, the left and right 
eigenvectors are identical and all eigenvalues are real. In this section we will not yet make 
this simplifying assumption. 



The potential of random frequencies can be treated as a perturbation of strength 
70" of the eigenvalue problem for the network Laplacian. Given the eigenvalues and 
orthonormal left and right eigenfunctions of L Eq. (|TTl) we are looking for the coefficients 
Eq^ of the expansion 

-Eo = Xo - laE^^^ - 7V^i?f - 0(7^=^) . (20) 

Again, we assume that the eigenvalue Aq = of the Laplacian is non-degenerate, so that the 
ground state is unique up to normalization. Accordingly, the synchronized state is unique up 
to a constant phase shift. Schrodinger perturbation theory, modified to allow for asymmetric 
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operators gives the expressions 

- E^') = (Po^ V,po) 



(21) 

Pi V,p,) (Pl Y,po] 



An — Ai- 

For the first two coefficients in the perturbation expansion Eq.Q of the synchronization 
frequency Q we find 



(22) 

Pj V.pfc) (Pl V,po) 



A, 



The ffist term is the weighted average of the frequencies with respect to the stationary 
probabihty distribution of the master equation P = L''"P with the transposed network 
Laplacian as matrix of transition rates. In the second expression we have used fi^^-* = —'^E^^ 
and Ao = 0. Equation Eq. (l22i) is a more compact form of equations Eqs. (IT^ - ffT^ combined. 



IV. Examples 



We can now study the change of the synchronization frequency with respect to oscillator 
heterogeneity and to the architecture of the coupling network. To find expressions for 
the expected first and second order response we will consider the ensemble of different 
realizations of random frequencies and an ensemble of random networks. Let us assume 
independent, identically distributed random frequencies with E [rjnrim] —^Vl?' = ^nm Then 
from Eq.(l22]) follows 



E p^^] = E, M 



(23) 



E 



Pl Vp„ p, 

Afr 



Here E^ is the expected value with respect to the frequency distribution and "^nw over 
the network ensemble. The vectors P^ and p^ in the second equation are left and right 



eigenvectors of the network Laplacian and 7^ the corresponding eigenvalues. The oper- 
ator is diagonal with the components of Pq on the diagonal. The first order frequency 
correction is independent of the network architecture while for uncorrelated frequencies the 
second order perturbation term is determined by the topology of the coupling network and 
the frequency heterogeneity a^. For symmetric coupling Anm = Amn the left eigenvector Pq 
is given by A^~^l and in the limit — ^ 00 the expected value of Vt^'^^ in Eq. (p^ can be 
written as 

E p'^] = 7 1 P(A) \ dX , (24) 
where p(A) is the Laplacian spectral density of the random network ensemble. 



her- 



32|. 



The integral Eq. ( 1241) has been studied in the different context of vibrational t 
modynamic stability for networks of linear springs, modelling complex molecules 
Whether the integral is finite depends on the spectral dimension d of the network, defined 
by the limit behavior of p(A) ~ A^^^ for A — 0, a suitable generalization of the Euclidean 



dimension for regular lattices to geometrically disordered structures [33|. For networks 
with spectral dimension d > 2 larger than two, the integral in Eq. fl24p is finite. In [32j] the 
authors study the case of a Sierpinski gasket which is a fractal graph for which the Laplacian 
spectrum can be calculated analytically and the spectral dimension is lower than two. In 



26l ] we study regular topologies for which the Fourier spectral decomposition is known and 
we also find the lower critical dimension d = 2. For ensembles of random graphs in general, 
it is a complicated task to find analytic expressions for the spectral density. Approximations 
of the spectral density of matrices associated with complex random networks, such as the 
Wigner semicircle law, are usually only available in the limit of dense networks, where the 
mean degree is much larger than one. 

A static scale-free random network model 



As an example we wil 
network model 



34 



use a recent result for the Laplacian of a static scale-free random 
35j. For this model the coupling strength Anm = ^mn between two 
oscillators is either zero or it is one with the probability kNwnWm, where the Wn ~ ■n,"^/^""^^ 
are normalized weights for the nodes n = 1 . . . N, and k is the mean degree of the network. 
The degree distribution follows a power law with exponent —a. In the thermodynamic limit 
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FIG. 1: (color online) Frequency shift in undirected, random networks of size N = 400 with 
Poissonian degree distribution and given mean degree k. The mean degree must be larger than 
kcr = 2 for an infinite random tree network. The frequency shift could be measured precisely 
by enforcing a zero mean frequency (D = 0. Random frequencies were drawn from a uniform 
distribution of standard deviation a = 10^^. From an ensemble of 10 realizations we show the mean 
network synchronization frequency divided by variance and nonisochronicity 7 (blue diamonds, 
Newton Method, 7 = 1.0), the predicted second order term {Q^'^^) from equation Eq. ()24p and the 
spectrum of the network Laplacian (red circles), the asymptotic line (dashed line) and the line 
{k — 2)~^ (solid line), which describes the actual behavior of the frequency shift even for small 
mean degrees close to fccr = 2. 



N CO and large A; 1 the spectral density of the Laplacian Eq.® is given [35| as 
P(A) = 



(a - 1)(-Ae)°-^ (-A)"" for A < Ac < , 
otherwise , 



(25) 



Ac = -k{a - 2){a - 1)-^ . 
Using this spectral density in equation Eq. (l21l) we obtain 



1 ._, (a-l)= 



E[n<'.] ^,/dAp(A)-^.r.l^. (26) 
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Erdos-Renyi model and random tree network limit 



The Erdos-Renyi random model 36| is recovered as a special case of the static scale-free 



random network model in the limit a oo 3J]. Then for ^ 1 we find E [f]*^^)] = 7^"^. 
One can study the limit of sparse uncorrelated random graphs by removing edges randomly 
without breaking the network in two components. The mean degree in a single component, 
undirected graph cannot be smaller than 2(A^ — 1)/N for a tree network. Every edge that 
is removed then breaks the connectivity, creates a new component and thus a new zero 
eigenvalue of the Laplacian. We expect therefore a divergence of the integral in Eq. flM|) 
for k ^ 2. We have tested our theory numerically for symmetrically connected random 
graphs with Poissonian uncorrelated degree distributions and = 400 phase oscillators 
with nonisochronicity 7 = 1.0. The frequencies were chosen randomly from a uniform 
distribution with var(c(j) = o"^ = 10^^. The synchronization frequency Q was determined 
on one hand by solving the algebraic equations Eq.Q with a Newton method, instead of 
integrating the KPEs Eq.Q, and on the other hand by using our perturbation approach 
and the complete eigenvalue spectrum of the network Laplacians. The results can be seen in 
Fig. ([T]). One can, indeed, see that the second order perturbation term diverges as {k — 2)~^ 
which is consistent with a powerlaw scaling E [fi*^^^] = 'jk~^ for larger mean degrees k. 

Application to network structure analysis 

In order to demonstrate possible applications of this perturbation theory to structural anal- 
ysis of an unknown coupling network let us now briefiy study what information can be 
gained from a measurement of linear and nonlinear responses to frequency changes of the 
oscillators. In 37|] the author presents a method to reconstruct a coupling network from 
measuring the linear response of the phase differences to linearly independent changes of 
the natural frequencies. This corresponds to using Eq. (JT3l) . The coupling network can be 



identified from the Green's function G of the network Laplacian and Eq. (fT3|) reads 

-d^^^ = Gt] . (27) 

If the phase differences are not accessible to direct measurement one can in principle also 
obtain the Green's function from the second order shift in the synchronization frequency. 
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For a symmetric coupling Eq. (122|) gives 

r](a;) = + 7 ^u;"^Glc; . (28) 

Let {a;^")} be a basis set of linear independent frequency detunings. Then the Green's 
function with respect to this basis can be determined from A^(A^ + l)/2 measurements of 
synchronization frequencies as 

Q + - {Q (a;(°)) + Q (a;^^))) = 27I u;'-"^^ Gu^^^ . (29) 

However, due to the number of measurements and the large time scales of a diffusion process 
the application is limited to very small networks, fast relaxation to the phase locked solution 
and high precision measurements. The analysis can be extended to nonidentical oscillators 
and asymmetric coupling. 



V. Discussion 



We have presented expressions for the first and second order perturbation terms of the 
synchronization frequency in complex networks of coupled Kuramoto phase oscillators with 
quenched frequency disorder. The two approaches in Sections II and III give equivalent 
results, but the second approach, based on a nonhnear approximation of the phase coupling 
function around zero, extends a well known treatment of the Kuramoto phase equations 



from continuous media to complex networks 17|, |26|, |27|. The results were given in terms 



of the eigenvalues and eigenvectors of the Laplacian matrix of the coupling network. In a 
sirigle component with mean degree of a undirected Erdos-Renyi random coupling network 



36( 1 and for oscillators with independent, identically distributed random frequencies Un of 
variance and nonisochronicity 7 the expected synchronization frequency was found to be 

E[n] = E[uj] + 7var(cj) {k - 2)-^ + 0(7V) . (30) 

While the expected synchronization frequency depends to the first order only on the natural 
frequencies in the system, the second order correction combines the nonlinearity 7 of the 
phase coupling function around zero, the variance of the frequencies and the mean degree 
of the coupling network in a simple way. 
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The explicit connection between synchronization frequency, natural frequencies and 
network structure in the Eq. fl22l) makes it in principle possible to infer information of 
either property from a measurement or the knowledge of the other properties. Network 
reconstruction by observing the linear response to a frequency detuning has already 



been proposed and successfully applied [37|. An analogous approach using frequency mea- 



surements instead of phase differences may be constructed based on the results of this paper. 

This work was supported by the German DFG through the project SfB555 and the 
Japanese JSPS. 

Appendix A 

Given the Kuramoto phase equations in synchronization 

N 

n = ar]n + '^ Anm 9{'&m " ^?n) , (31) 

and a phase locked solution of the KPEs for identical oscillators 

N 

^^'^ = Y.^nrng{^t^-^T), (32) 
m=l 

we want to derive expressions for the coefficients in the expansion of the synchronization 
frequency Vt in powers of the frequency heterogeneity a. Let us start by implicitly defining 
notations for the involved perturbation terms and phase differences 



n = + , (33) 

1=1 

oo 

= ^i°)+^n = (34) 

1=1 

^mn = '&m — = ^mn + V^mn = ^mn + "^mn 5 (35) 



1=1 



oo 



9mn = 9 i^^Z) , 9^ = d^9mn ■ (37) 
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Note that here we do not assume = const, g{0) = or gmn = g'^^^O). It has been 
pointed out, that even for identical oscillators the homogeneous solution may not be the 
only synchronized solution of the Kuramoto phase equations Iss!]. In certain coupling 
topologies and for large nonisochronicity the completely synchronized solution can coexist 



with a dominating chaotic attractor of drifting phases [23|]. If the network is homogeneous 
and sufficiently well connected, however, the stable solution of complete synchronization is 
typical. Therefore we assume = g{0) = and grnn — g^'^\0) in the main text of this 
paper. 



The coefficients (finn yield the recursion relation 



rmn 



•dmn for j = 1 , 

Et'i^Mn'"'"'^ for j</, (38) 
otherwise . 



Inserting Eq. llMI) into Eq. (l3T]) we find 



N 



Q = ar]n + ^ Anm gmn , (39) 

m=l 

N oo ^ 

= +(Tr]n + Y, ^nm -gi 



I Jmnrmn ' 

m=l j=l 



N oo 



, - , r^mn 5 

?! — ' 

m=l j=l l=j 



-l,l~k) 



oo I N ^ l-l 

1=1 j=l m=l •'' k=l 

In the second line we have inserted the unperturbed solution Eq. fl32l) and in the third line 
we used the expansion Eq. (l36|) of the powers of (fmn- Since the leading order of ipmn in a is 
one, the jth power has a leading term of order j. In the last line the recursion relation fl55]) 
was used. If we now collect the nonlinear terms (j > 1) in a vector b^'^ we can write down 
this result in a more compact form 

- = E ^'^^'^ = E f E ^nmgLn^^l + , (40) 

1=1 1=1 \m=l j 
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or in vector form 

n«i= (L^(') + b«) , (41) 

where 1 is a constant vector of unit entries, i?*-'-* is the vector of perturbative corrections 
to the phases and b*^') is a vector which only depends on perturbation terms of order lower 
than I. The matrix L is the Jacobian 

N 

Lnm — ^nm 9mn ^nm ^ ^ ^nl 9ln ' (^^) 

1=1 

or the Laplacian if g'^^ is a constant. The vectors b*^'-' are given by 

6« = Vn , (43) 

6(2) = \^ A -a" t^W^ (U) 

m=l 

&(3) = (a" , 1 /// ^{i)3\ /^gv, 

m=l ^ ^ 

and in general for / > 1 

IN l-l 
j=2 m=l k=l 

Equation Eq. fj4T|) can be solved iteratively for each perturbation order. Let us consider a 
complete, orthonormal set of left and right eigenvectors and of the Jacobian with 

Lpfc = AfcPfc , L^Pfc = A^Pfc , 

(47) 

N-l 

Pfc Pfc' = ^kk' , PfcPfc = I 5 



fc=0 



and in particular 

po = 1 , r Po = 1 . (48) 

We can now define the projectors 

Po = PoPS, Qo = I-Po. (49) 
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The operation Pqx projects to a constant vector where all entries are equal to the weighted 
average {x)pg and Qo removes this average from the components of a vector. Applying these 
projectors to Equation Eq. lHTj) we obtain 

Q(0 = pt b(0 ^ (50) 

= L^(') + Qob(') . (51) 
The last equation is solved for i?'-'-' up to an arbitrary global phase shift by 



pt b(0 

^^'^ = -EM-^P^- (52) 
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